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The  design  of  the  flow  channels  in  PEM  fuel  cells  directly  impacts  the  transport  of  reactant  gases  to  the 
electrodes  and  affects  cell  performance.  This  paper  presents  results  from  a  study  to  optimize  the  geometry 
of  the  flow  channels  in  a  PEM  fuel  cell.  The  optimization  process  implements  a  genetic  algorithm  to 
rapidly  converge  on  the  channel  geometry  that  provides  the  highest  net  power  output  from  the  cell.  In 
addition,  this  work  implements  a  method  for  the  automatic  generation  of  parameterized  channel  domains 
that  are  evaluated  for  performance  using  a  commercial  computational  fluid  dynamics  package  from 
ANSYS.  The  software  package  includes  GAMBIT  as  the  solid  modeling  and  meshing  software,  the  solver 
FLUENT,  and  a  PEMFC  Add-on  Module  capable  of  modeling  the  relevant  physical  and  electrochemical 
mechanisms  that  describe  PEM  fuel  cell  operation.  The  result  of  the  optimization  process  is  a  set  of 
optimal  channel  geometry  values  for  the  single-serpentine  channel  configuration.  The  performance  of 
the  optimal  geometry  is  contrasted  with  a  sub-optimal  one  by  comparing  contour  plots  of  current  density, 
oxygen  and  hydrogen  concentration.  In  addition,  the  role  of  convective  bypass  in  bringing  fresh  reactant 
to  the  catalyst  layer  is  examined  in  detail.  The  convergence  to  the  optimal  geometry  is  confirmed  by  a 
bracketing  study  which  compares  the  performance  of  the  best  individual  to  those  of  its  neighbors  with 
adjacent  parameter  values. 

©  201 1  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

A  polymer  electrolyte  membrane  fuel  cell  (PEMFC)  is  a  device 
that  converts  the  chemical  energy  of  hydrogen  directly  into  elec¬ 
trical  energy.  The  main  components  of  a  PEMFC  are  assembled 
in  a  layered  configuration  as  shown  in  Fig.  1.  The  central  layer 
is  the  proton  exchange  membrane  (PEM)  which  is  coated  with 
catalyst  layers  on  either  side  that  serve  as  the  sites  for  the  oxi¬ 
dation  (anode)  and  reduction  (cathode)  reactions.  On  either  side  of 
the  catalyst-coated  membrane  are  gas  diffusion  layers  (GDL)  com¬ 
prised  of  porous  carbon  paper  that  are  used  to  deliver  the  reactant 
gases  (hydrogen  and  oxygen)  to  the  electrodes  from  the  supply 
channels  in  the  bipolar  plates  which  are  the  outermost  layers  of 
the  assembly.  The  flow  channels  simultaneously  remove  product 
water  from  the  fuel  cell.  The  bipolar  plates  also  collect  current 
from  the  fuel  cell  and  serve  as  the  anode  and  cathode  terminals. 
They  are  typically  made  of  graphite,  composites,  or  stainless  steel 
for  considerations  of  cost,  structural  rigidity,  thermal  conductiv¬ 
ity,  and  electrical  conductivity.  Depending  on  the  material  of  the 
bipolar  plate,  the  flow  channels  are  machined,  molded  or  stamped. 
Ribs  located  between  adjacent  channels  are  responsible  for  mak¬ 
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ing  electrical  contact  with  the  GDL,  which  is  itself  in  direct  contact 
with  the  electrode.  The  contact  region  between  the  ribs  and  the 
GDL  is  known  as  the  land  area.  The  GDL  is  typically  0.2-0.3  mm 
thick  and  serves  the  purpose  of  evenly  distributing  the  reactants 
to  the  catalyst  layers  while  also  creating  electrical  contact,  and 
cushioning  the  membrane  electrode  assembly  (MEA)  against  the 
clamping  pressure  under  the  lands. 

The  anode  and  cathode  half-cell  reactions  are  shown  in  Fig.  1. 
Although  the  ideal  voltage  for  this  type  of  fuel  cell  is  1 .23  V,  various 
loss  mechanisms  reduce  the  practical  operating  voltage  to  about 
0.6  V.  The  main  loss  mechanisms  include  activation  loss,  ohmic  loss, 
and  mass  transport  loss.  Activation  loss  is  a  result  of  slow  reaction 
kinetics,  especially  at  the  cathode.  Ohmic  loss  is  due  primarily  to 
the  resistance  to  proton  flow  in  the  membrane,  as  well  as  the  elec¬ 
trical  resistance  of  the  bipolar  plates  and  the  various  interconnects 
within  the  cell.  Mass  transport  loss  is  typically  noticed  at  high  cur¬ 
rent  draws  due  to  fuel  starvation  and  possible  flooding.  Running 
the  fuel  cell  at  higher  stoichiometries  can  alleviate  mass  transport 
losses,  but  at  the  cost  of  higher  parasitic  pumping  losses. 

The  focus  of  this  research  is  on  the  geometric  design  of  the 
bipolar  plates  of  the  cell.  As  stated  earlier,  flow  channels  within 
these  plates  deliver  the  reactants  to  the  catalyst  layers,  while  also 
removing  the  products.  The  design  of  the  flow  channels  directly 
impacts  the  transport  of  reactant  gases  to  the  anode  and  cath¬ 
ode  catalyst  layers  and  strongly  affects  cell  performance.  A  large 
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Fig.  1-  Schematic  cross  section  of  a  PEMFC. 


focus  on  channel  design  can  be  found  in  the  literature.  The  basic 
channel  designs  include  serpentine,  parallel,  and  interdigitated 
flow  fields.  Each  basic  design  has  its  own  strengths  and  weak¬ 
nesses.  Parallel  channels  experience  the  smallest  pumping  losses, 
but  since  the  pressure  drop  across  the  length  of  the  channels  is 
small,  they  are  prone  to  blockage  by  liquid  water  and  regions 
of  the  cell  become  incapable  of  producing  current.  Interdigitated 
channels  are  more  effective  at  water  removal  and  have  enhanced 
convective  bypass  into  the  GDL  [1],  but  experience  higher  pres¬ 
sure  drops  across  the  cell  since  the  inlet  and  outlet  are  not  directly 
connected  by  channels.  The  serpentine  design  is  quite  common 
(see  Shimpalee  et  al.  [2])  because  it  is  more  effective  at  expelling 
liquid  water  thereby  reducing  the  problem  of  channel  flooding. 
In  addition,  the  pressure  difference  between  adjacent  channels 
enhances  convective  bypass  into  the  GDL  and  improves  reactant 
delivery  to  the  catalyst  layer  [1  ].  However,  the  serpentine  channel 
experiences  higher  pumping  losses  due  to  its  high  overall  channel 
length. 

Each  type  of  channel  design  is  characterized  by  several  param¬ 
eters  that  control  gas  transport  and  impact  the  cell’s  performance. 
For  channels  of  rectangular  cross  section,  these  parameters  include 
the  channel  width,  channel  height,  land  width,  and  channel  length. 
It  is  important  to  determine  how  these  geometric  parameters  influ¬ 
ence  fuel  cell  performance  in  order  to  design  efficient  fuel  cells. 
The  selection  of  values  for  these  parameters  has  been  the  subject 
of  many  research  investigations  and  has  been  explored  experi¬ 
mentally,  numerically,  and  analytically.  For  example,  a  study  on 
land-to-channel  width  ratio  was  done  by  Yan  et  al.  [4].  Tapering  the 
width  continuously  between  the  inlet  and  outlet  was  investigated 
by  Yan  et  al.  [5].  Xu  and  Zhao  [6]  developed  a  convection-enhanced 
serpentine  flow  field  (CESFF)  by  forcing  a  single  channel  to  double 
back  on  itself.  All  of  these  examples  identify  the  optimal  param¬ 
eters  for  a  cell  of  given  dimensions  and  tend  to  generalize  that 
these  parameters  can  be  used  to  design  fuel  cells  of  larger  dimen¬ 
sions. 

Performance  improvements  obtained  by  designing  cells  that 
make  highly  effective  use  of  a  given  nominal  area  can  lower 
the  cost  of  fuel  cells  which  will  promote  their  commercializa¬ 
tion.  The  goal  of  the  current  work  is  to  devise  a  method  that 
can  rapidly  optimize  the  geometric  design  of  a  flow  field  that 
maximizes  the  net  power  output  for  an  active  area  of  given 
dimensions.  Our  optimization  method  employs  a  genetic  algo¬ 
rithm  to  efficiently  determine  the  optimal  channel  geometry 
within  a  relatively  large  search  space.  The  effectiveness  of  a  spe¬ 
cific  fuel  cell  is  evaluated  by  a  fitness  function  whose  value  is 
the  net  power  per  unit  area,  or  net  power  density.  A  set  of  in- 
house  scripts  was  developed  to  combine  the  strengths  of  finite 
element  analysis  (FEA)  with  a  parameterized  solid  model  draw¬ 
ing  script  and  a  genetic  algorithm  optimization  method.  The 
channel  designs  are  considered  for  a  specific  nominal  area  and 


for  a  specific  set  of  operating  conditions.  The  principal  contri¬ 
bution  of  this  work  is  the  automated  optimization  method  to 
efficiently  determine  the  best  channel  design  for  a  PEMFC  of  a  given 
size. 


2.  PEMFC  model 

In  order  to  minimize  the  costs  associated  with  the  design  and 
optimization  of  an  individual  cell,  it  is  valuable  to  attain  predictive 
capability  regarding  the  cell’s  performance  potential  before  actual 
construction  and  testing.  It  is  therefore  beneficial  to  include  numer¬ 
ical  modeling  in  the  PEMFC  design  process.  Many  works  have  been 
dedicated  to  the  construction  of  numerical  models  that  can  accu¬ 
rately  predict  PEMFC  performance.  These  models  have  increased 
in  complexity  over  time,  evolving  from  two-dimensional  models 
by  Chen  et  al.  [7],  to  relatively  complex  three-dimensional  mod¬ 
els.  Such  complex  models  have  been  created  by  different  groups 
including  Le  and  Zhou  [8 ]  who  included  multiphase  flow,  Shimpalee 
et  al.  [2],  and  Fan  et  al.  [9].  These  researchers  have  constructed  their 
own  in-house  CFD  models  or  used  commercially  available  codes. 
For  our  analysis  we  use  a  commercially  available  software  package 
(FLUENT  6.3.26). 


2.1.  Model  assumptions 

The  ANSYS  software  package  FLUENT  and  the  FLUENT  PEM  fuel 
cell  module  is  used  in  this  research  to  compile  the  appropriate 
user-defined  functions  for  a  PEMFC.  This  model  does  not  take  into 
account  areas  of  importance  such  as  material  degradation,  the  pres¬ 
ence  of  impurities  in  the  reactant  supply  streams,  and  non-ideal 
operating  conditions  such  as  variances  in  control  system  responses 
and  the  non-isothermal  nature  of  the  environment  which  the  fuel 
cell  is  bound  to  experience  during  operation.  Each  of  those  factors 
would  lead  to  reduced  fuel  cell  performance.  Despite  these  limita¬ 
tions  on  the  absolute  accuracy  of  the  predicted  performance,  it  is 
sufficient  to  obtain  the  relative  performance  of  the  various  cell  con¬ 
figurations  for  the  purposes  of  this  study.  Likewise,  product  water 
is  assumed  to  remain  exclusively  in  vapor  state  in  this  model.  This 
limitation  is  also  acceptable  for  our  purpose  as  the  chosen  operat¬ 
ing  condition  is  sufficiently  removed  from  the  flooding  regime.  The 
PEMFC  model  used  in  this  research  assumes  that  the  flow  is  laminar, 
the  fluid  is  incompressible,  the  inlet  gases  follow  the  ideal  gas  law, 
the  gas  diffusion  layers,  catalyst  layers,  and  membrane  layer  are 
isotropic  materials,  the  flow  is  single  phase,  and  that  steady-state 
conditions  exist.  The  experimental  results  obtained  by  Spernjak 
et  al.  [10]  in  the  Fuel  Cell  Research  Laboratory  at  the  University 
of  Delaware  will  be  used  to  verify  the  PEMFC  model’s  results  in  a 
later  section. 
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2.2.  Governing  equations 

Reactant  transport  is  strongly  influenced  by  the  geometry  of  the 
reactant  flow  channels  in  a  PEMFC.  The  important  governing  equa¬ 
tions  are  described  below,  beginning  with  the  conservation  of  mass 
equation: 

^  +  V.(/ov)  =  Sm  (1) 


The  conservation  of  species  is  required  to  calculate  the  mass  bal¬ 
ance  for  each  of  the  involved  reactants  in  this  model  and  is  defined 
below: 

+  v  •  (pvyj)  =  V  ■  (pDjVyj)  +  Sj  (8) 

where  Sj  refers  to  the  species  source  term. 

The  effective  mass  diffusivity  is  calculated  using  the  expression 
by  [11,12]. 


Eq.  (1 )  is  valid  for  mass  balance  in  the  channels  as  well  as  the  GDL 
where  both  diffusion  and  convection  are  important.  Sm  is  the  mass 
source  term  and  for  all  zones  in  this  model  Sm  ~  0. 

Next,  the  momentum  equation  is  used  to  solve  for  the  fluid 
velocities  in  the  channels  and  the  GDL  and  the  species’  partial  pres¬ 
sures: 

^  +  V-(/ow)  =  -Vp  +  V-r  +  SM  (2) 

The  momentum  source  term,  5m,  is  zero  in  all  zones  except  for  the 
gas  diffusion  layers  and  the  catalyst  layers.  In  these  zones,  the  Darcy 
equation  is  used  to  define  5M. 

Sm  =  -^£V  (3) 

Here,  fi  is  the  viscosity  of  the  fluid,  and  I<  and  s  are  the  permeability 
and  porosity,  respectively,  of  the  particular  zone. 

Temperature  is  modeled  through  the  energy  equation  whose 
general  form  is  as  follows: 


(9) 


The  species  source  terms,  Sj,  are  zero  except  for  the  fluid  zones  and 
are  described  by  the  next  three  equations: 


_  Mw,  h2 
^h2  -  2^Ka 

Mw,02  n 

•J02  — - TE — Kcat 


2  F 


Sh2o  =  -- 


Mi 


w,H20  . 

—  1 
2  F 


(10) 

(11) 

(12) 


The  electrochemistry  modeling  begins  by  calculating  the  rates 
of  hydrogen  oxidation  and  oxygen  reduction  in  the  anode  and  cath¬ 
ode  catalyst  layers,  respectively.  The  following  are  the  two  potential 
equations  solved  in  the  PEM  model.  Eq.  (13)  is  solved  for  the  elec¬ 
tron  transport  through  the  solid  conducting  regions  including  the 
bipolar  plates,  gas  diffusion  layers  and  catalyst  layers.  Eq.  (14)  is 
solved  for  the  proton  conductivity  through  the  membrane. 


ft{eph  +  (\-e)pshs]  +  V  (pvh)  =  V-  [  keffVT-J2hjJj  ]  +sh 


(4) 


where  h  is  the  enthalpy  of  the  gas  mixture,/  is  the  diffusive  flux,  the 
subscript  j  represents  the  chemical  species  of  the  bracketed  quan¬ 
tity,  and  the  subscript  s  represents  the  solid  phase  of  the  bracketed 
quantity.  Zones  with  solid  phases  include  the  GDL,  CL,  membrane, 
and  the  bipolar  plate.  The  first  two  terms  on  the  right  hand  side  of 
Eq.  (4)  are  the  conduction  and  species  diffusion  terms,  respectively. 
Sft  is  the  volumetric  source  term  and  is  defined  for  all  zones  of  the 
cell  by: 


Sh  —  I  Rohm  +  h 


\ironrt\f\in  I  nRn 


(5) 


where  I2R0hm  is  the  ohmic  heating  term,  hreaction  is  the  heat  of  for¬ 
mation  of  water  term,  q R  is  the  electric  work  term  {q  and  R  are 
the  activation  loss  and  volumetric  transfer  current,  respectively), 
and  /iphase  is  the  latent  heat  of  water  term.  The  effective  thermal 
conductivity  is  calculated  as: 


keff  =  skf  +  (1  -e)ks  (6) 

where  the  subscripts /and  s  represent  the  fluid  and  solid  compo¬ 
nent,  respectively,  of  a  specific  zone. 

The  heterogeneous  reactions  that  take  place  on  the  catalyst  sur¬ 
faces  are  balanced  by  their  rate  of  production  as  described  below: 


pDy 

5 


M, 


( yj.surf  -  yj.centfr  = 


(7) 


where  Dj  is  the  effective  mass  diffusivity,  r  is  the  surface  to  volume 
ratio  of  the  catalyst  layer,  8  is  the  average  distance  between  reaction 
surfaces  and  cell  center,  y  is  the  mass  fraction  of  species,  and  the 
subscripts  j,surf  and  /cent  refer  to  the  surface  and  the  center  of 
the  cell,  respectively.  MWJ  is  the  molecular  weight  of  the  species,  n 
represents  the  number  of  electrons  transferred  per  molecule  and  F 
is  the  Faraday  number. 


V-(osV0s)  +  Ks  =  0  (13) 

V  •  (<XmV0m)  +  Rm  =  0  (14) 

In  the  preceding  equations,  crs  is  the  electrical  conductivity  of  the 
respective  solid  zone,  om  is  the  membrane  electrical  conductivity 
calculated  by  Eq.  (15)  [12],  0  is  the  electric  potential,  and  R  is  the 
volumetric  transfer  current. 

am  =  /3s( 0.514A.  -  o.26)"e1268(1/303-1/r)  (15) 


Here,  A  and  co  are  water  content  and  the  membrane  conductiv¬ 
ity  constant,  respectively.  For  the  anode  side,  the  source  terms  are 
defined  as  Rs  =  -Ran  and  Rm  =  +Ran-  On  the  cathode  side,  the  source 
terms  are  defined  as  Rs  =  +Rcat  and  Rm  =  - RCat •  The  source  terms  for 
the  unnamed  zones  are  assumed  to  be  equal  to  zero. 

The  terms  Ran  and  Rcat  are  defined  according  to  the 
Butler-Volmer  equations  as  shown  below. 


Van 


^^^2  J  ^gUanFhanlRT  _  g—oicatFhanlRT^ 


Ycat 


Rcat  =  j'ft  J  (-e'Ja"F’lmtlRT  +  e-«cntF(?cnt/RT) 


(16) 


(17) 


The  local  activation  losses  are  solved  as  follows: 


qan  —  <Ps  ~  0m  (18) 

qcat  =  0.s  ~  0m  —  (19) 

The  activation  loss,  q,  is  calculated  as  the  difference  between  the 
solid  and  membrane  potentials  as  shown  above. 


2.3.  PEMFC  material  properties  and  boundary  conditions 

This  section  summarizes  the  boundary  conditions  and  material 
properties  that  are  maintained  constant  in  the  current  simulations. 
The  inlet  gases  are  set  at  a  constant  mass  flow  rate  of  humidified 
hydrogen  (anode  side)  and  air  (cathode  side).  The  mass  flow  rates 
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Table  1 

Material  properties. 


Property 

Value 

Property 

Value 

Cp\col,gdl,cl 

871 J  kg-1  K"1 

&an,cat 

2 

Cp-,mem 

2000 J  kg-1  K"1 

P 

1 

D° 

j 

3  x  10“5  m2  s_1 

Van, cat 

1 

hgdi 

0.25  mm 

Yp 

1 

hd 

12.5  |jim 

Yt 

1.5 

hmem 

25  fxm 

£cl,mem 

0.5 

jref 

Jan 

5  x  109  Am-3 

Pcol,cl,gdl 

2719  kg  m~3 

fef 

Jcat 

4  x  109  Am-3 

Pmem 

1980  kg  m~3 

kp;col,gdl,cl 

8  Wm1  K"1 

&cl,gdl 

5000  x  107  ohm-1  m_1 

kgdl.cl 

1  x  10“12  m2 

&COI 

1.0  x  106  ohm-1  rrr1 

Mm 

llOOkgkmoU1 

&  mem 

1.0  x  10-16  ohm-1  nr1 

Po 

101.325  kPa 

CO 

1 

r 

4  x  109  rrr1 

T0 

300  K 

Xj,ref 

1 

Table  2 

Boundary  conditions. 

Location 

Boundary  condition  type 

Value 

Cathode  and  anode  inlets 

Mass  flow  rate 

[Calculated] 

Anode  and  cathode  outlets 

Pressure 

101325  Pa 

All  exterior  surfaces 

Temperature 

353 1< 

Cathode  terminal  surface 

Voltage 

0.6  V 

Anode  terminal  surface 

Voltage 

0.0  V 

are  calculated  for  a  stoichiometry  of  A,  =  2  at  1 A  cm-1.  Using  this 
condition,  appropriate  values  for  the  mass  fluxes  are  calculated 
based  on  the  area  of  a  given  cell.  The  outlet  of  each  channel  is  open 
to  the  atmosphere  and  set  to  1  atm.  The  cathode  terminal  wall  is  set 
at  a  constant  0.6  V  and  the  anode  terminal  wall  is  grounded  at  0  V, 
and  both  terminal  wall  surfaces  are  set  to  80  °C.  Tables  1  and  2  sum¬ 
marize  the  material  properties  and  boundary  conditions  employed 
in  this  PEMFC  model. 

2.4.  Meshing  strategy  and  convergence  considerations 

Building  a  three-dimensional  solid  model  with  an  appropriately 
constructed  mesh  is  an  important  concern  when  performing  any 
type  of  numerical  modeling  and  simulation  routine.  The  mesh  used 
here  corresponds  to  the  best  practice  recommendations  of  the  FLU¬ 
ENT  PEMFC  Add-on  Module.  Each  zone  of  the  cell  requires  a  specific 
set  of  meshing  rules.  The  GDL  thickness  is  divided  into  1 0  mesh  lay¬ 
ers,  the  catalyst  layer  and  membrane  thicknesses  are  divided  into 
four  mesh  layers,  and  the  channel  height  is  divided  into  10  lay¬ 
ers.  The  bends  of  each  serpentine  pass  have  the  highest  velocity 
and  pressure  gradients,  so  it  is  important  to  keep  the  mesh  den¬ 
sity  high  in  these  zones,  but  the  mesh  density  is  allowed  to  coarsen 
by  a  symmetric  1.1  meshing  ratio  along  each  pass.  The  result  is 
that  the  channel  ends  are  defined  by  cube  shaped  mesh  elements, 
while  the  channel  lengths  are  defined  by  slightly  elongated  hexag¬ 
onal  elements.  An  important  feature  of  the  meshing  process  used 
by  the  Gmesh  script  ensures  that  there  are  no  skewed  elements. 
All  meshed  elements  are  hexagonal  with  90°  angles  between  each 
face.  An  example  mesh  is  shown  in  Fig.  2. 

A  mesh  refinement  study  was  conducted  in  order  to  show  that 
the  solution  was  grid  independent.  In  particular,  the  mesh  density 
was  increased  at  the  channel  ends  where  the  velocity  gradients  are 
highest.  The  final  mesh  density  selected  for  the  simulations  showed 
negligible  changes  in  the  solution  with  further  refinement. 

The  solver  is  based  on  the  finite-volume  method  and  employs 
a  pressure-based  solver  with  double-precision.  The  convergence 
of  the  model  depends  on  several  constraints.  First,  all  the  residual 
quantities  including  continuity,  three  velocity  components,  energy, 
hydrogen  species,  oxygen  species,  water  species,  electrical  poten¬ 


Fig.  2.  Example  of  mesh  density  in  the  single-serpentine  fuel  cell  domain  (top  bipo¬ 
lar  plate  removed  in  order  to  reveal  the  channel  domain). 

tial,  protonic  potential,  and  water  content  must  attain  a  tolerance  at 
or  below  1  x  10~4.  The  convergence  criteria  also  require  the  average 
y-current  density  on  the  cathode  and  the  pressure  at  the  cathode 
inlet  to  converge  to  below  1  x  10-6.  The  y-current  density  is  the  cur¬ 
rent  density  normal  to  the  surface  of  the  bipolar  plate.  These  two 
quantities  are  required  to  determine  the  gross  power  density  of  the 
cell,  and  the  parasitic  power  consumption.  These  values  are  needed 
in  the  optimization  strategy  and  will  be  discussed  in  Section  6. 

Our  optimization  approach  requires  multiple  simulations  with 
varying  geometrical  characteristics  for  the  flow  channels.  Varying 
the  channel  geometry  requires  the  construction  of  a  unique  model 
and  mesh  for  each  simulation.  Therefore,  an  efficient  and  method¬ 
ical  way  of  constructing  and  meshing  a  wide  variety  of  PEMFCs 
is  required  for  this  study.  Drawing  and  meshing  each  unique  cell 
would  require  a  prohibitively  large  amount  of  time  and  therefore, 
the  process  was  automated  using  a  GAMBIT  journal  file,  which  will 
be  referred  to  hereon  as  the  Gscript  file.  GAMBIT  [13]  is  a  com¬ 
mercially  available  software  package  that  was  used  in  this  study  to 
create  the  domain  of  a  PEMFC,  divide  the  model  into  a  discretized 
mesh,  and  identify  the  necessary  boundary  conditions  and  zones 
without  manual  input.  GAMBIT  outputs  a  mesh  file  that  can  be  read 
by  FLUENT  [14]  to  perform  a  simulation  at  the  desired  operating 
conditions. 

For  the  case  of  a  single-serpentine  flow  field  design,  four  param¬ 
eters  are  needed  to  fully  describe  the  design  for  a  fixed  area: 
channel  width  (CW),  channel  height  (CH),  land  width  (LW),  and 
channel  length  (CL).  The  Gscript  file  receives  these  quantities  as 
input  parameters  and  appropriately  constructs  the  solid  model. 
However,  in  this  study,  CH  is  set  to  a  constant  value  of  1.0  mm. 

3.  PEMFC  model  validation  and  verification 

3.1.  Comparison  with  experimental  data 

In  order  to  validate  the  model,  a  comparison  to  an  experimen¬ 
tal  result  is  provided  in  this  section.  The  model  was  constructed  to 
replicate  the  cell  employed  by  Spernjak  et  al.  [  1 0]  who  used  1 0  cm2 
PEMFC  bipolar  plates  made  of  graphite  with  a  single-serpentine 
channel  configuration.  The  channels  were  0.8  mm  wide,  1  mm  high, 
and  the  land  width  was  0.8  mm.  Humidified  air  and  hydrogen 
were  supplied  at  0.18  and  0.076  slpm,  respectively.  The  cell  was 
set  to  operate  at  a  constant  temperature  of  80  °C.  The  OCV  for 
this  verification  run  was  set  to  0.778  V  to  match  the  correspond¬ 
ing  experimental  value.  The  simulation  was  initiated  by  setting 
the  voltage  to  OCV  followed  by  stepwise  0.5  V  decrements  corre¬ 
sponding  to  the  experimental  data.  Fig.  3  shows  a  good  agreement 
between  the  experimental  and  numerical  results.  The  small  differ¬ 
ences  between  the  model  and  the  experimental  data  are  due  to  the 
assumptions  made  in  the  model.  However,  the  overall  degree  of 
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agreement  demonstrates  that  the  model  is  a  good  approximation 
over  the  prevailing  range  of  operating  conditions. 

In  addition,  global  mass  conservation  was  also  checked  to  verify 
the  model.  This  check  was  accomplished  by  examining  simulation 
data  for  one  point  on  the  polarization  curve  corresponding  to  the 
operating  voltage  of  0.624  V.  FLUENT  calculates  the  amount  of  oxy¬ 
gen  consumed  by  the  cell  and  the  average  current  flux  magnitude 
in  the  y-direction  at  the  cathode  terminal  wall.  The  total  release  of 
electrons  can  be  directly  related  to  the  oxygen  consumption  rate  as 
given  below. 

,  mo24F  (6.629  x  10-4)(4)(9.6485  x  104)  ,  nnc  A 

I=~M~  = - 3T9988 - =  7.995  A  (20) 

Here,  fn  is  the  oxygen  consumption  rate  (g  s-1 ),  M  is  the  molecular 
mass  of  oxygen,  and  F  is  the  Faraday  constant.  Next,  the  current 
density  over  the  cathode  terminal  surface  was  integrated  to  obtain 
a  total  current  of  7.993  A,  which  agrees  very  well  with  the  value  in 
Eq.  (20)  and  confirms  that  mass  is  globally  conserved. 


Fig.  4.  Comparison  of  numerical  and  analytical  PEMFC  pressure  distribution  along 
the  center  pass. 

provides  further  evidence  that  the  model  is  capable  of  accurate 
predictions. 

4.  Optimization  strategy 

The  genetic  algorithm  (GA)  is  an  optimization  technique  that 
employs  the  fundamental  principles  that  drive  natural  evolution 
and  was  originally  developed  by  Holland  according  to  Sivanandam 
and  Deepa  [15].  The  strength  of  the  GA  is  its  ability  to  efficiently 
locate  optima  within  relatively  large  search  spaces.  In  order  for  a 
GA  to  be  implemented,  a  system  must  first  be  reduced  to  the  basic 
parameters  that  are  needed  to  describe  it.  Once  these  parameters 
are  identified,  the  GA  can  go  to  work  on  determining  the  optimum 
values  for  these  parameters.  The  GA  used  in  this  study  is  the  com¬ 
mercially  available  Genetic  Algorithm  and  Direct  Search  Toolbox 
available  from  MatLab. 

4.1.  Genetic  algorithm  theory 


3.2.  Comparison  with  analytical  results 

Although  the  experimental  verification  shows  that  the  model  is 
capable  of  predicting  overall  performance  values  in  the  form  of  cur¬ 
rent  density  output  versus  voltage,  it  is  necessary  to  further  verify 
detailed  results  internal  to  the  cell.  Feser  et  al.  [3]  formulated  an 
analytical  expression  that  predicts  the  pressure  distribution  along 
the  channels  of  the  cell: 

In  this  equation,  A PceU  is  the  total  pressure  drop  across  the  entire 
cell,  Nc  is  the  number  of  passes  in  the  cell,  CL  is  the  length  of  a 
single  pass,  and  m  is  a  dimensionless  quantity  defined  as  m2  = 
4(L2A“1}(hGDL/LWX/qp//<c/ian).  In  this  expression,  Ac  is  the  prod¬ 
uct  of  the  channel  width  and  channel  height,  kip  is  the  in-plane 
permeability  of  the  GDL,  and  kchan  is  the  permeability  of  the  chan¬ 
nel.  The  channel  permeability  is  approximated  as  a  function  of  the 
height-to-width  ratio  of  the  channel  according  to  Feser  et  al.  [3]. 

Fig.  4  provides  a  comparison  between  the  pressure  distribution 
along  the  center  pass  predicted  by  the  model  with  the  analytical 
values  provided  by  Eq.  (21 ).  The  data  in  Fig.  4  have  been  normalized 
by  the  length  (CL)  of  a  single  pass  and  by  the  maximum  pres¬ 
sure  along  the  pass  ( Pmax )•  The  agreement  between  these  results 


The  description  of  the  GA  in  this  research  employs  terms  typi¬ 
cally  used  when  discussing  genetics.  An  individual  is  defined  by  the 
genes  that  make  up  its  chromosome.  In  this  work,  the  chromosome 
is  responsible  for  describing  a  unique  PEMFC  individual  (IND).  Each 
individual  has  a  specific  number  of  genes.  For  a  three-parameter 
optimization,  there  are  three  genes,  one  for  each  parameter.  The 
parameters  in  this  research  include  the  channel  width,  land  width, 
and  channel  length;  each  parameter  can  assume  values  within 
an  allowable  range  that  is  set  by  selecting  an  upper  and  lower 
bound. 

The  GA  proceeds  by  spawning  new  generations,  such  that  indi¬ 
viduals  in  the  current  generation  are  created  by  examining  and 
selecting  genes  from  the  individuals  in  the  previous  generation. 
The  performance  of  an  individual  is  determined  by  applying  an 
appropriately  formulated  fitness  function  in  the  GA.  The  fitness 
function  is  responsible  for  minimizing  or  maximizing  a  quantity 
of  interest.  In  the  case  of  the  PEMFC,  the  fitness  function  is  defined 
as  the  net  power  density.  The  GA  includes  the  crossover  process 
between  successive  generations  such  that  the  genes  of  the  best  indi¬ 
viduals  are  crossed  with  those  of  other  individuals  which  showed 
good  performance.  The  new  individuals  are  comprised  of  differ¬ 
ent  combinations  of  the  genes  in  the  previous  generation.  The  GA 
also  includes  mutations  between  successive  generations  to  help 
the  optimization  process  by  including  values  that  were  not  present 
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Table  3 

Genetic  algorithm  settings. 


GA  setting 

Chosen  setting 

Crossover  fraction 

0.8 

Stall  time  limit 

Infinite 

TolFun 

1  x  10“6 

TolCon 

1  x  10“6 

Creation  function 

Uniform 

Fitness  scaling  function 

Rank 

Selection  function 

Stochastic  uniform 

during  earlier  generations.  Random  mutations  help  to  dislodge  a 
solution  that  may  have  settled  on  a  local  optimum  and  propel  it 
towards  the  desired  global  optimum.  The  mutation  rate  is  set  so  as 
to  not  create  completely  random  individuals,  but  rather  to  include 
a  limited  amount  of  randomness  between  successive  generations. 
While  the  method  does  not  guarantee  that  the  global  optimum  will 
always  be  found,  it  is  generally  accepted  that  the  final  solution  will 
be  close  to  the  global  optimum  as  long  as  successive  generations 
are  incapable  of  producing  better  individuals.  Not  finding  the  exact 
global  optimum  is  a  trade  off  for  saving  computational  time  and 
power. 

The  GA  finds  the  optimum  solution  with  far  fewer  simulations 
than  an  exhaustive  brute  force  parametric  study  which  evaluates 
the  performance  of  each  of  a  vast  array  of  possible  individuals.  The 
MatLab  GA  and  Direct  Search  Toolbox  are  capable  of  executing  such 
a  scheme  along  with  many  other  options  that  are  necessary  to  con¬ 
trol  a  GA  optimization.  This  toolbox  is  fully  integrated  with  the 
automated  process  that  optimizes  a  PEMFC  in  this  research.  The 
following  sections  will  outline  the  options  that  guide  the  GA  and 
the  integration  of  each  of  the  software  tools  that  are  required  to 
automate  the  optimization  and  solution  processes. 

4.2.  MATLAB  GA  Toolbox  set  parameters 

Table  3  summarizes  the  GA  parameters  that  are  used  through¬ 
out  this  study.  The  elite  count  is  the  number  of  individuals  that  are 
preserved  from  one  generation  to  the  next.  These  elite  individu¬ 
als  are  the  best  performers  of  a  given  generation  and  are  used  to 
crossover  and  create  the  rest  of  the  population  in  a  given  genera¬ 
tion.  The  crossover  rate  can  range  from  zero  to  one  and  represents 
the  fraction  of  the  next  generation  that  is  comprised  of  crossover 
offspring.  The  remaining  offspring  are  comprised  of  mutated  indi¬ 
viduals.  The  stall  time  limit  is  set  to  infinity  in  order  to  eliminate 
the  possibility  of  terminating  the  PEMFC  model  before  complet¬ 
ing  a  simulation.  The  tolerance  values  are  set  to  terminate  the  GA 
process  if  successive  generations  do  not  create  new  individuals 
whose  performance  improvement  exceeds  the  set  tolerance  value. 
This  setting  would  suggest  that  the  optimum  was  reached  and  the 
individuals  converged  to  the  best  solution  within  the  search  space. 
The  fitness  scaling  function  puts  the  individuals  of  a  given  popula¬ 
tion  in  rank  order  instead  of  ranking  them  by  their  fitness  values. 
The  best  individual  is  ranked  number  1,  and  so  on.  The  ranking  of 
individuals  removes  the  effects  of  the  spread  between  raw  fitness 
scores. 

4.3.  Communication  script 

Now  that  the  GA  settings  have  been  fully  defined  and  the  opti¬ 
mization  process  has  been  described,  the  process  of  evaluating  an 
individual’s  fitness  needs  to  be  discussed.  This  study  represents  the 
first  instance  in  which  GA  optimization  has  been  coupled  with  an 
automated  scheme  capable  of  constructing  unique  PEMFC  mod¬ 
els  and  evaluating  their  performance.  Although  previous  studies 
have  investigated  PEMFC  performance  parametrically,  none  have 
been  aimed  at  finding  the  optimal  channel  characteristics  of  a 


customizable  area.  The  following  script  is  capable  of  automati¬ 
cally  generating  a  solid  model  for  a  fuel  cell  of  a  chosen  design 
and  size.  The  channel  geometry  that  produces  optimum  perfor¬ 
mance  for  a  cell  of  a  certain  size  need  not  be  the  same  as  that 
for  a  larger  cell.  The  benefit  of  our  script  is  the  scalability  of 
the  process;  an  optimum  geometry  can  be  found  for  a  cell  of 
any  size,  automatically  and  efficiently.  This  research  identifies  the 
optimum  values  for  a  three-parameter  single  serpentine  configu¬ 
ration. 

Several  scripts  were  written  to  control  and  integrate  all  of  the 
software  packages  necessary  to  run  the  optimization.  Fig.  5  summa¬ 
rizes  the  flow  and  connectivity  of  the  various  scripts.  The  user  sets 
the  number  of  parameters  to  be  optimized  as  well  as  the  upper  and 
lower  bounds  for  each  parameter.  Once  the  previously  discussed  GA 
options  are  set  and  the  search  space  is  chosen,  the  GA  optimization 
is  ready  to  begin. 

At  this  point,  the  GA  Toolbox  creates  random  sets  of  param¬ 
eter  values  for  each  of  the  individuals  in  the  initial  population. 
This  initial  population  is  referred  to  as  generation  (GEN)  zero.  The 
GA  evaluates  the  fitness  of  each  individual  by  calling  the  fitness 
function  (FF).  The  fitness  function  is  the  script  that  was  created  to 
control  the  optimization  and  is  therefore  responsible  for  the  inte¬ 
gration  of  all  of  the  software  packages  and  extracts  the  appropriate 
data  to  calculate  the  fitness  of  an  individual.  The  script  calculates 
the  fitness  of  each  individual  of  the  generation,  one  at  a  time.  The 
fitness  function  receives  as  input  the  values  for  each  of  the  genes 
that  describe  an  individual.  Each  of  the  values  is  written  to  a  GAM¬ 
BIT  journal  file.  The  GAMBIT  journal  file,  Gmesh,  is  the  third  script 
in  this  process.  It  is  designed  to  draw  a  solid  model  according 
to  the  input  parameter  values.  It  then  populates  the  model  with 
an  appropriately  scaled  mesh  and  applies  the  boundary  condition 
types  and  zone  types  to  the  PEMFC  domain.  Once  the  model  is  fully 
constructed,  the  mesh  is  written  to  file. 

The  FF  script  next  calls  the  FLUENT  journal  script,  Fjou.  The  Fjou 
script  is  responsible  for  loading  the  PEMFC  module,  applying  all 
material  types  to  each  of  the  zones,  applying  values  to  each  of  the 
required  boundary  conditions,  initializing  the  solver,  setting  con¬ 
vergence  criteria,  converging  the  solution,  and  saving  the  desired 
solution  values  to  a  file.  The  Gmesh  script  takes  the  raw  data  from 
the  Fjou  output  files  and  uses  it  to  calculate  the  fitness  of  each  indi¬ 
vidual.  This  final  value  is  returned  to  FF  and  sent  back  to  the  GA 
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control  process.  This  process  is  repeated  for  each  of  the  individuals 
that  the  GA  requests  for  each  generation. 

Specific  strengths  of  the  scripts  are:  the  Gmesh  script  calculates 
the  number  of  serpentine  passes  that  are  required  to  fill  a  given 
area;  the  Fjou  script  can  handle  any  Gmesh  script  with  appropri¬ 
ately  named  zones;  and  the  operating  conditions  are  easy  to  access 
if  alterations  are  desired. 

The  GA  terminates  when  the  specified  number  of  generations  is 
reached,  or  when  the  convergence  criteria  for  successive  genera¬ 
tions  are  met.  The  GA  then  outputs  the  best  individual  of  the  entire 
population  along  with  its  associated  parameter  values.  This  individ¬ 
ual  represents  an  optimum  configuration  within  the  search  space. 
The  individual  can  then  be  compared  to  other  individuals  in  the 
population  to  extract  the  characteristics  that  strongly  influenced 
their  performance.  The  information  gathered  from  the  population 
as  a  whole  can  be  used  to  design  and  initialize  future  GA  studies  in 
order  to  decrease  convergence  time  and  increase  performance. 


5.  Results  and  discussion 

This  study  uses  a  population  of  eight  individuals  for  each  gen¬ 
eration,  and  a  limit  of  12  new  generations  beyond  the  initially 
generated  randomized  population  as  summarized  in  Table  3.  Each 
new  generation  selects  two  elite  individuals  from  the  previous 
generation.  These  two  elite  individuals  are  used  to  create  four 
additional  crossover  offspring.  The  remaining  two  individuals  of 
each  generation  are  created  by  the  mutation  function.  The  entire 
population  consists  of  104  individuals;  however,  since  the  elite 
individuals  pass  unchanged  between  generations,  the  number  of 
unique  individuals  that  need  to  be  evaluated  is  equal  to  80.  The 
fitness  function  script  prevents  the  reevaluation  of  individuals  by 
checking  a  continually  updated  database  for  repeat  individuals  dur¬ 
ing  the  optimization. 

The  search  space  for  the  design  of  this  PEMFC  is  defined  by  pre¬ 
scribing  ranges  of  values  for  the  three  independent  parameters: 
channel  length  (CL),  channel  width  (CW),  and  land  width  (LW). 
The  channel  length  is  defined  as  the  length  of  a  single  pass  of  the 
serpentine  channel,  as  opposed  to  the  overall  length  of  the  entire 
serpentine  channel  from  inlet  to  outlet.  The  three  parameters  are 
depicted  schematically  in  Fig.  6. 

The  active  area  of  the  fuel  cell  extends  beyond  the  outer  edges 
of  the  channel  in  the  xz-plane  by  an  amount  equal  to  one-half 
of  the  specified  land  width  as  shown  in  Fig.  6.  The  nominal  area 
(A)  of  the  cell  is  set  at  a  target  value  of  16  cm2.  The  range  of  val¬ 
ues  for  each  of  the  parameters  is  as  follows:  4cm  <  CL <  14cm; 
0.5  mm  <  CW  <1.5  mm;  and  1 .0  mm  <  LW  <3.5  mm.  The  channel 
height  is  set  to  a  constant  value  of  1.0  mm.  The  Gmesh  script  uses 
the  three  parameter  values  to  calculate  the  number  of  passes  that 
are  required  to  fill  the  nominal  area  with  the  given  channel  width, 
land  width,  and  channel  length.  The  number  of  passes  is  rounded 
to  the  nearest  odd  integer  in  order  to  create  a  cell  with  the  desired 
channel  configuration  in  which  the  inlet  and  outlet  ports  are  on 
opposite  sides  of  the  cell.  This  condition  maintains  the  same  inlet 
and  outlet  configuration  for  all  individuals.  The  PEMFC  is  designed 
with  counter-current  flow  channels  such  that  the  inlet  and  outlet 
ports  of  the  anode  are  reversed  for  the  cathode.  The  anode  and  cath¬ 
ode  channels  overlay  exactly  on  top  of  each  other  when  projected 
on  to  the  xz-plane. 

The  gross  power  density  of  a  cell  is  the  product  of  its  current  den¬ 
sity,  i,  and  its  voltage,  V.  In  practice,  a  fuel  cell  system  incorporates 
the  balance-of-plant  consisting  of  several  components  whose  input 
power  is  drawn  from  the  fuel  cell  itself.  The  cumulative  balance-of- 
plant  power  requirement  represents  the  cell’s  parasitic  power  loss. 
Of  particular  relevance  to  this  study,  the  power  requirement  of  the 
motor  and  compressor  assembly  for  the  air  supply  system  is  largely 


influenced  by  the  channel  geometry.  Therefore,  it  is  appropriate  to 
include  the  parasitic  loss  associated  with  driving  air  through  the 
cathode  channels  as  a  factor  in  the  fitness  function.  Each  channel 
configuration  within  the  search  space  imposes  a  power  require¬ 
ment  that  is  proportional  to  the  total  pressure  drop  from  inlet  to 
outlet  and  the  mass  flow  rate  through  the  channel.  A  small  chan¬ 
nel  cross-sectional  area  and  long  overall  channel  length  results  in 
a  much  higher  pressure  drop  than  a  channel  with  a  large  cross- 
sectional  area  and  short  overall  channel  length.  Larminie  and  Dicks 
[16]  provide  an  expression  for  the  parasitic  power  requirement  of 
a  cell: 


In  this  equation,  the  combined  motor  and  compressor  efficiency 
r\  =  0.7,  and  the  heat  capacity  ratio  of  air  y  =  1 .4.  Tamb  is  the  ambient 
temperature,  pamb  is  the  ambient  pressure,  pin  is  the  pressure  at  the 
inlet  of  the  PEMFC,  and  m  is  the  mass  flux  of  air.  The  inlet  pressure 
is  calculated  by  FLUENT  and  is  part  of  the  data  output  of  the  Fjou 
script.  The  mass  flow  rate  of  air  is  calculated  using  an  expression 
from  Larminie  and  Dicks  [16]: 

mair  =  3.57  x  10-7  •  XiA  (23) 

Here,  i  is  the  current  density,  and  A  is  the  nominal  area  of  the  cell. 
Similarly,  the  hydrogen  flow  rate  is  calculated  as: 

mH2  =  1.05  x  10“8  XiA  (24) 

Both  reactant  gas  mass  flow  rates  are  set  to  provide  a  stoichiom¬ 
etry  of  A  =  2  at  /  =  lAcm-2,  which  is  a  common  condition  used 
when  evaluating  PEMFC  performance.  With  all  of  the  required 
quantities  defined,  the  net  power  is  now  calculated  in  Eq.  (25). 
The  parasitic  power  term  does  not  include  any  loss  due  to  the 
supply  of  hydrogen  to  the  anode  side  because  the  hydrogen  is  sup- 
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Table  4 

Best  and  worst  individual  comparison. 


Quantity  or  parameter 

Best 

Worst 

Channel  length 

4.64  cm 

4.0  cm 

Channel  width 

0.762  mm 

1.11  mm 

Land  width 

0.934  mm 

3.50  mm 

Net  power  density 

0.954  W  cm"2 

0.474  W  cm-2 

plied  from  a  pressurized  vessel,  and  therefore  requires  no  external 
compression. 

Pnet  —  Pcell  ~  Ppar  (25) 

The  power  of  the  cell,  Pceu,  is  calculated  by  the  PEMFC  module  and 
is  the  product  of  the  cell’s  current  density  integrated  over  the  nom¬ 
inal  area  and  the  cell  voltage.  The  desired  fitness  function  is  finally 
defined  as  the  net  power  per  unit  area. 


The  negative  sign  is  included  in  the  fitness  function  because  the 
MatLab  GA  technique  searches  for  the  minimum  in  the  fitness  func¬ 
tion,  whereas  our  goal  is  to  maximize  the  net  power  density.  Once 
the  simulation  is  complete,  the  absolute  value  of  the  fitness  func¬ 
tion  represents  the  net  power  density  of  each  individual. 

5. I .  Comparison  of  best  and  worst  individuals 

This  GA  optimization  strategy  resulted  in  a  total  of  104  individ¬ 
uals,  80  of  which  were  unique.  These  individuals  covered  a  wide 
range  of  parameter  values  within  the  limits  set  in  the  GA.  The  per¬ 
formance  corresponds  to  the  fitness  function  defined  in  Eq.  (26). 
Because  the  fitness  function  is  the  net  power  density  as  opposed  to 
the  gross  power  density,  the  individuals  are  ranked  more  accurately 
on  how  they  would  perform  during  actual  operation. 

The  focus  of  this  section  is  on  the  comparison  between  the  best 
and  worst  individuals  in  the  population  whose  parameter  and  out¬ 
put  values  are  summarized  in  Table  4.  This  comparison  aids  in  the 
explanation  of  the  physical  reasons  for  why  two  individuals  within 
a  reasonable  search  space  produce  very  different  performance  val¬ 
ues.  The  net  power  density  of  the  worst  individual  is  50.7%  lower 
than  that  of  the  best  individual. 

The  investigation  of  these  individuals  begins  by  examining  the 
current  density  distributions  on  the  surface  of  the  cathode  catalyst 
layer  of  each  individual.  Fig.  7  shows  the  current  density  magnitude 


on  the  cathode  GDL/CL  interface  where  the  oxygen  reduction  reac¬ 
tion  takes  place.  Areas  of  relatively  high  current  density  indicate 
a  more  vigorous  reaction  and  good  oxygen  transport.  The  figure 
shows  that  the  worst  individual  is  unable  to  utilize  the  catalyst 
sites  under  the  land  regions  of  the  cell,  whereas  the  best  individual 
makes  much  better  use  of  the  areas  under  the  land.  In  the  worst 
individual,  the  majority  of  the  current  is  produced  directly  under 
the  channel  and  is  collected  by  the  ribs  of  the  bipolar  plate.  The 
reason  for  the  poor  utilization  of  active  area  under  the  lands  in  the 
worst  individual  is  the  high  land  width  of  3.5  mm.  Reactant  gas  is 
unable  to  diffuse  laterally  far  under  the  land  and  results  in  mass 
transport  loss  leading  to  poor  current  density.  It  is  also  evident  that 
the  inlet  and  outlet  regions  do  not  produce  much  current  density 
in  either  of  the  individuals.  This  is  due  to  the  counter-current  flow 
configuration  used  in  the  cell.  As  a  result,  by  the  time  the  reactant 
flow  reaches  the  outlet,  reactant  depletion  lowers  its  partial  pres¬ 
sure,  which  starves  the  region  of  either  oxygen  or  hydrogen  and 
therefore  limits  the  reaction. 

Next,  we  investigate  the  oxygen  consumption  within  the  cell. 
When  a  PEMFC  cathode  is  fed  with  air  instead  of  pure  oxygen, 
the  reaction  is  more  likely  to  be  limited  by  the  oxygen  content 
of  the  air  between  the  inlet  and  the  outlet  of  the  cell.  To  under¬ 
stand  a  given  cell’s  performance,  it  is  important  to  examine  the 
local  consumption  of  oxygen  and  correlate  it  to  the  overall  perfor¬ 
mance  of  the  cell.  Fig.  8  compares  the  best  and  worst  individuals 
in  terms  of  the  oxygen  concentration  in  the  cathode  channel 
midplane. 

Fig.  8  shows  that  for  the  best  individual  the  oxygen  is  almost 
completely  depleted  by  the  time  the  flow  reaches  the  outlet.  This 
observation  confirms  that  oxygen  consumption  is  high  for  the  best 
individual  resulting  is  high  overall  current  density.  On  the  other 
hand,  the  worst  individual’s  outlet  oxygen  concentration  is  still 
quite  high,  about  half  of  the  inlet  concentration.  The  main  reason 
for  this  is  the  reduced  oxygen  consumption  due  to  the  channel’s 
inability  to  deliver  reactants  to  the  catalyst  layer  underneath  the 
wide  land  regions  via  convective  bypass. 

Convective  bypass  is  the  primary  mechanism  that  delivers  reac¬ 
tant  gases  to  the  catalyst  layer  under  the  lands.  Convective  bypass 
is  proportional  to  the  pressure  difference  between  two  adjacent 
channels  in  a  cell.  The  largest  pressure  difference  occurs  between 
the  ends  of  two  adjacent  channels  across  the  root  of  the  rib  (across 
points  A  and  B  in  Fig.  6).  If  sufficiently  large,  this  pressure  differ¬ 
ence  can  convectively  drive  reactant  gas  into  the  GDL  and  the  CL 
under  the  lands,  while  also  driving  out  product  water  from  under 
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Fig.  7.  Contour  plots  of  current  density  (Am-2)  at  the  cathode  GDL/CL  interface  for  the  best  (left)  and  worst  (right)  individuals. 
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Fig.  8.  Contour  plots  of  oxygen  concentration  (kmol  m-3)  along  the  cathode  channel  for  the  best  (left)  and  worst  (right)  individuals. 


the  lands.  According  to  Darcy’s  Law,  convective  bypass  through  the 
GDL  is  inversely  proportional  to  the  width  of  the  land  region.  Hence, 
a  large  land  width  can  severely  inhibit  convective  bypass  even  in  the 
presence  of  a  significant  pressure  differential.  The  large  land  width 
of  3.5  mm  for  the  worst  case  individual  is  responsible  for  extended 
dead  zones  under  the  lands  and  poor  overall  performance. 

In  addition,  the  reaction  can  proceed  only  if  there  is  an  ade¬ 
quate  proton  flux  through  the  membrane  from  the  anode  side  to 
complete  the  reaction  on  the  cathode  side.  If  the  hydrogen  con¬ 
centration  is  depleted  such  as  near  the  hydrogen  outlet,  then  not 
enough  hydrogen  is  convectively  forced  into  the  GDL  and  the  reac¬ 
tion  rate  and  current  production  are  weak.  The  occurrence  of  such 
dead  zones  should  be  minimized  as  much  as  possible.  Fig.  9  shows 
that  there  is  still  a  significant  concentration  of  hydrogen  at  the 
outlet  of  the  worst  performer,  while  the  hydrogen  is  almost  com¬ 
pletely  depleted  by  the  time  the  flow  reaches  the  outlet  of  the  best 
individual. 

The  contour  plot  of  oxygen  concentration  in  three  layers  (top, 
middle  and  bottom)  of  the  cathode  GDL  in  Fig.  1 0  shows  that  there  is 
oxygen  under  both  the  channel  and  the  land  for  the  best  performer. 
This  figure  also  shows  that  oxygen  is  driven  towards  the  catalyst 
layer  by  the  concentration  gradient  of  oxygen  in  the  y-direction. 
Oxygen  diffuses  through  the  thickness  of  the  GDL  parallel  to  the 


y-axis  and  is  consumed  directly  below  the  channels.  In  addition, 
oxygen  is  also  being  delivered  to  the  catalyst  under  the  land  areas. 

The  velocity  vectors  in  Fig.  1 1  provide  evidence  that  air  is  driven 
by  convection  under  the  ribs  of  adjacent  passes  for  the  best  per¬ 
former.  Therefore,  oxygen  is  delivered  effectively  to  the  catalyst 
layer  even  below  the  land  areas.  This  is  a  desirable  result  because 
all  of  the  catalyst  area  is  able  to  participate  in  the  reaction.  In 
contrast,  the  velocities  under  the  cathode  lands  in  the  case  of 
the  worst  individual  are  not  very  high,  implying  that  diffusion  is 
the  dominant  transport  mechanism  for  oxygen.  Diffusion  is  not  as 
effective  in  transporting  reactant  gas  under  the  lands,  and  there¬ 
fore,  although  the  oxygen  concentration  is  high  under  the  channels 
for  the  worst  individual,  it  experiences  low  oxygen  concentration 
under  the  lands.  Contrary  to  this  result,  the  best  performer  shows 
evidence  of  both  diffusive  and  convective  transport. 

The  absolute  pressure  distribution  within  the  channels  and  the 
GDL  provides  insight  into  the  velocity  distributions  found  in  the 
GDL.  Fig.  12  presents  the  pressure  distribution  along  the  midplane 
of  the  cathode  GDL  for  the  best  and  worst  performers.  As  expected, 
the  pressure  gradients  are  highest  in  the  vicinity  of  the  channel 
roots  at  each  pass.  In  the  best  individual,  the  pressure  gradient  is 
very  effective  at  driving  convective  bypass  under  the  land  areas  of 
the  cell.  The  worst  individual  suffers  because  the  pressure  gradient 
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Fig.  9.  Contour  plots  of  hydrogen  concentration  (kmol  m-3 )  along  the  anode  channel  for  the  best  (left)  and  worst  (right)  individuals. 
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Fig.  10.  Oxygen  concentration  (kmol  m-3)  at  the  top,  middle,  and  bottom  of  the  cathode  GDL  for  the  best  (top)  and  worst  (bottom)  performers.  The  bottom  of  the  GDL  is  in 
contact  with  the  catalyst  layer. 


is  insufficient  to  drive  convective  bypass  across  the  wide  lands.  The 
overall  pressure  difference  across  the  worst  case  cell  is  small  owing 
to  the  small  number  of  passes,  and  hence  the  small  total  channel 
length.  The  pressure  distribution  is  an  important  indicator  of  the 
cell’s  performance.  The  pressure  distribution  is  mainly  determined 
by  the  configuration  of  the  gas  flow  channels.  This  result  further 
validates  the  usefulness  of  this  study  to  determine  the  optimal 
channel  configurations  for  a  specified  cell  area. 


5.2.  Comparison  to  neighbors 

The  data  presented  thus  far  have  shown  that,  of  a  given  set  of 
individuals,  there  is  clearly  one  that  performs  the  best.  The  pre¬ 
ceding  discussion  has  also  helped  to  explain  the  mechanisms  that 
contribute  to  its  good  performance.  Keeping  in  mind  that  the  set 
of  individuals  that  was  actually  evaluated  by  the  GA  is  just  a  small 
subset  of  all  the  possible  individuals  the  search  space  encompasses, 


Fig.  11.  Velocity  vectors  (ms-1)  along  the  midplane  of  the  cathode  GDL  for  the  best  (top)  and  worst  (bottom)  individuals. 
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Fig.  12.  Absolute  pressure  (Pa)  along  the  midplane  of  the  cathode  GDL  for  the  best  (left)  and  worst  (right)  performers. 


it  is  necessary  to  provide  more  evidence  that  the  “best”  individ¬ 
ual  is  indeed  the  best  of  all  possible  individuals.  The  only  way  to 
prove  this  conclusively  is  to  perform  an  exhaustive  study  of  all  of 
the  combinations  of  the  parameter  values.  This  type  of  brute-force 
study  is  prohibitively  expensive,  and  the  GA  technique  was  cho¬ 
sen  in  the  first  place  to  precisely  avoid  such  an  exercise.  A  simpler 
approach  is  to  examine  neighboring  individuals  with  parameter 
values  adjacent  to  those  of  the  best  individual  and  thereby  confirm 
that  the  best  individual  is  superior  to  its  neighbors.  Ideally,  the  GA 
should  have  ignored  such  inferior  individuals  and  converged  to  the 
best  individual.  Accordingly,  a  “bracketing”  study  was  performed 
in  which  each  parameter  value  was  varied  from  the  optimum  value 
by  ±5%  and  ±10%  while  keeping  the  other  parameter  values  fixed. 
The  bracketing  study  resulted  in  12  new  individuals  for  whom  the 
corresponding  fitness  values  were  evaluated.  The  performance  of 
the  new  individuals  is  compared  with  the  best  performer  in  Fig.  13. 

Fig.  13  shows  that  while  all  of  the  new  individuals  performed 
well,  none  was  able  to  out-perform  the  best  individual.  This  plot  has 
two  interesting  characteristics.  First,  the  GA  did  indeed  correctly 


Fig.  13.  Result  of  a  bracketing  study  that  compares  the  performance  of  neighbors 
with  that  of  the  best  individual. 


identify  the  best  individual  within  this  parametric  “bracketing” 
study.  Second,  the  best  performer  occupies  the  peak  of  a  relatively 
flat  performance  region.  The  neighbors  surrounding  the  best  indi¬ 
vidual  show  a  performance  drop  of  no  more  than  2%  of  the  optimal 
configuration.  This  observation  suggests  that  while  performing  the 
study  on  this  size  of  a  PEMFC,  it  is  important  to  isolate  channel 
parameters  that  are  near  the  optimum  value,  and  that  being  close 
to  the  best  performer  is  also  a  very  good  result. 

5.3.  Parasitic  pumping  losses 

Upon  inspecting  Fig.  14,  it  becomes  obvious  why  the  fitness 
function  was  chosen  to  include  the  parasitic  power  losses  asso¬ 
ciated  with  driving  air  through  the  cell.  The  parasitic  power  loss 
does  not  vary  monotonically  with  the  gross  power  output  of  the 
cell;  i.e.,  a  ranked  order  of  best  performers  by  net  power  would  not 
be  identical  to  a  similar  list  ranked  by  gross  power.  This  is  because 
for  a  given  cell  area,  the  channel  configuration  can  take  on  various 
channel  geometries.  The  narrower  the  channel  cross-sectional  area 


Fig.  14.  Comparison  of  gross  and  net  power  density  between  the  best  individual 
and  its  neighbors. 
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and  the  longer  the  overall  channel  length,  the  higher  the  pumping 
requirements  of  the  cell.  Therefore,  it  is  necessary  to  account  for  the 
parasitic  power  losses  to  accurately  gauge  PEMFC  performance. 

6.  Conclusions 

This  research  has  focused  on  the  optimization  of  a  three- 
parameter  single-serpentine  channel  configuration  using  the 
genetic  algorithm  method.  The  best  possible  individual,  identi¬ 
fied  as  the  candidate  with  the  highest  net  power  density,  was 
automatically  determined  by  the  GA  and  reasons  for  its  superior 
performance  were  discussed.  The  GA  was  proven  to  effectively 
search  a  space  for  the  best  individual  by  means  of  a  bracketing  study 
which  demonstrated  that  the  best  performer  is  indeed  superior  to 
all  its  immediate  neighbors  in  parameter  space. 

Through  this  research,  a  powerful  method  of  incorporating  the 
strengths  of  several  commercially  available  software  packages  has 
been  developed.  The  unique  contribution  of  this  work  is  the  method 
to  efficiently  automate  the  optimization  of  a  PEMFC  for  a  specified 
application.  This  method  has  developed  an  interface  to  automate 
communication  of  inputs  and  outputs  between  GAMBIT,  FLUENT, 
and  the  MatLab  GA  Toolbox.  Each  software  package  has  a  unique 
scripting  language  capable  of  producing  the  data  flow  required 
by  the  optimization  process.  The  communication  is  accomplished 
through  the  combination  of  several  scripts,  written  in-house, 
that  utilize  these  coding  languages  as  well  as  the  command 
line. 

GA  optimization  has  proven  to  be  an  efficient  tool  for  searching 
a  relatively  large  search  space  for  an  optimal  value.  The  algo¬ 
rithm  was  successful  in  correctly  identifying  the  best  individual 
even  though  the  search  space  associated  with  this  work  yielded  a 
relatively  flat  region  for  the  value  of  the  fitness  function  in  its  imme¬ 
diate  neighborhood.  Therefore,  the  optimization  can  be  trusted  to 
find  an  optimum  relatively  close  to  the  global  search  space  optimal 
value. 

Some  limitations  of  this  project  include  the  following  important 
considerations.  First,  the  PEMFC  model  employed  here  assumed 
single-phase  flow.  Multiphase  flow  was  not  considered  in  this 
PEMFC  model  in  order  to  reduce  model  complexity  and  compu¬ 
tational  time.  Second,  the  best  individuals  tended  to  consume  all 
of  the  available  reactants  implying  some  degree  of  reactant  starva¬ 
tion  near  the  outlets.  These  cells  could  have  benefited  from  a  higher 
inlet  mass  flow  rate  such  that  the  reactant  utilization  was  not  quite 
so  high,  as  long  as  the  increased  pumping  power  did  not  cause  the 
net  power  density  of  the  cell  to  drop. 

While  the  search  space  explored  by  this  algorithm  thus  far  is 
useful  for  lab-scale  fuel  cell  evaluation  and  comparison,  the  real 
value  is  in  its  ability  to  find  an  optimal  PEMFC  design  configura¬ 


tion  of  any  size.  By  successfully  demonstrating  the  optimization  of 
laboratory-scale  PEMFCs,  the  method  can  confidently  be  applied  to 
fuel  cells  of  a  larger  target  area.  It  is  expected  that  in  larger  cells,  the 
physical  mechanisms  that  govern  PEMFC  performance  would  yield 
different  optimal  parameter  values.  The  method  developed  here  is 
capable  of  including  as  many  parameters  as  the  user  desires,  includ¬ 
ing  material  properties.  Such  additional  parameters  may  include, 
but  are  not  limited  to:  channel  height,  GDL  thickness,  GDL  porosity, 
GDL  permeability,  and  channel  taper. 

A  significant  feature  of  the  method  developed  here  is  that  other 
types  of  flow  fields  can  also  be  evaluated.  Configurations  such  as 
parallel,  interdigitated,  and  multiple  serpentine  channels  can  be 
incorporated  using  the  Gmesh  script  developed  for  this  study.  As 
long  as  the  same  zones  and  boundary  condition  types  are  prescribed 
by  the  Gmesh  script,  the  Fjou  script  is  capable  of  evaluating  the  fit¬ 
ness  of  any  PEMFC  and  its  desired  channel  configuration.  However, 
it  is  recommended  that  mesh  refinement  studies  are  conducted 
on  new  channel  configurations  to  ensure  proper  convergence  and 
solution  accuracy. 
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